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Although geometrical frustration transcends scale, it has primarily been evoked in the micro 
and mesoscopic realm to characterize such phases as spin-ice liquids and glasses and to explain 
the behavior of such materials as multiferroics, high temperature superconductors, colloids and 
copolymers. Here we introduce a system of macroscopic ferromagnetic rotors arranged in a planar 
lattice capable of out-of-plane movement that exhibit the characteristic honeycomb spin ice rules 
studied and seen so far only in its mesoscopic manifestation. We find that a polarized initial state 
of this system settles into the honeycomb spin ice phase with relaxation on multiple time scales. We 
explain this relaxation process using a minimal classical mechanical model which includes Coulombic 
interactions between magnetic charges located at the ends of the magnets and viscous dissipation at 
the hinges. Our study shows how macroscopic frustration arises in a purely classical setting that is 
amenable to experiment, easy manipulation, theory and computation, and shows phenomena that 
are not visible in their microscopic counterparts. 



Frustration in physical systems commonly arises be- 
cause geometrical or topological constraints prevent 
global energy minima from being realized. Although not 
limited to microscopic phenomena, it is commonly seen 
in compounds with spins forming lattices with a trian- 
gular motif In such systems, frustration may lead to 
the existence of ice selection rules Q which have been 
observed in a variety of materials where spins form net- 
works such as the corner-sharing tetrahedra, known as 
the Pyrochlore lattice 0-111 ; leading to monopole-like ex- 
citations Q and other exotic phases of matter Q. Even 
though, artificial spin ices [8l-ll0| have shown that frustra- 
tion can be mimicked by classical magnets, these systems 
do not account quantitatively for the effects of inertia, 
dissipation Ill4l3j , dilution and geometrical disorder be- 
cause of the mesoscopic scale and fast dynamics of the 
domain walls (~ 10 ns) that hinder the understanding of 
collective dynamics processes. Here we aim to circum- 
vent this situation by introducing a new macroscopic re- 
alization of a frustrated magnetic system created using 
single out-of-plane rotational degree of freedom magnetic 
rotors, arranged in a kagomc lattice, a pattern of corner- 
sharing triangular plaquettes that dynamically evolves 
into a spin-ice phase after a magnetic quench. The ice 
phase is reached due to the delicate interplay between in- 
ertia, friction and Coulomb-like interactions between the 
macroscopic magnetic rods. Our prototypical frustrated 
system has a few advantages for research in frustrated 
magnetic systems associated with the ability to (i) tune 
the interactions through changes in distance and/or ori- 
entation between magnets and (ii) examine the lattice 
relaxation dynamics by direct visualization at a single 
particle level. 

A minimal macroscopic realization of local frustration 
can be seen easily in a 120° star configuration using three 
ferromagnetic rods with their hinges on a plane (Fig. [TJt) . 
The rods have length L = 2a = 1.9 x 10~ 2 m, diameter 
d = 1.5 X 1(T 3 m, mass M = 0.28 x 10~ 3 Kg and satu- 



ration magnetization M s — 1.2 x 10 6 A m _1 . By design 
the only allowed motions for the rotors are rotations in 
the polar direction a. The hinges supporting the rods 
were placed at the sites of a kagome lattice with lattice 
constant I = \/3(a + A) where A is the shortest dis- 
tance between the tips and the nearest vertex center and 
A/L ~ 0.2 (Fig. OJ,), so that when in the x~y plane, the 
magnets realize the bonds of a honeycomb lattice. The 
magnetization of a rotor i is defined as the vector mi join- 
ing its N to its S pole, thus mi is the coarse-grained spin 
variable for each magnet. When all three magnets are 
close to each other, the lowest energy configuration con- 
sists of one pole being different from the others, leading 
to a frustrated state consisting of permutations of NNS or 
SSN (S=south pole, N = north pole) that correspond to 
the honeycomb spin ice rules 0, HH . With this unit-cell 
plaquette, we prepare a polarized lattice of n = 352 of 
these magnetic rotors, with an unavoidable geometrical 
disorder in the azimuthal orientation of the rotors, 6, due 
to lattice imperfections 59 max ~ 2°; this follows a Gaus- 
sian distribution with mean 58 = 1.2°. We oriented the 
S poles of all rotors out of the plane by applying a strong 
magnetic field along the z direction B z = 3.2 x 10~ 3 T 
(Supplementary Information S4). At t — the field was 
switched off, to allow for the lattice to relax, a procedure 
that was repeated several times. After about 2 seconds, 
all the rotors had reached equilibrium configurations very 
near the x — y plane (the non-planarity out of the x — y 
plane 5a ~ 10° on average) and in the honeycomb spin 
ice manifold. Figure shows a picture of the lattice 
where all the rods fulfill the ice rules. The experimental 
distribution of vertices is shown in the (red) bars of Fig- 
ure [Tfc. We find all vertices falling into the six low energy 
(spin ice) configurations while high energy states (type 1 
and 2) are absent. 

This macroscopic spin ice consists of elemental rotor 
units that constitute a frustrated triad which we char- 
acterize at a static and a dynamic level (Supplementary 
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Information SI, S2 and S3). This allows us to use a dipo- 
lar dumbbell approach to the magnets 0] , determine the 
charge q = irM s d 2 /4 ~ 2.03 A m, at each pole, find the 
damping time scale for an isolated rod t-q ~ 1 s and 
examine how Coulomb interactions and geometrical dis- 
order in 9 and A control the orientations of the rods rel- 
ative to each other. On a collective level, the relaxation 
of the lattice from the z polarized state to the spin ice 
manifold may be characterized in terms of the correlation 
between nearest neighbor spins a and /3, with S a Sp = 1 
when m Q ■ is positive, S a S/3 = —1 otherwise. From 
high speed movies (400 fps), we extracted the full time 
trajectory ai(t) of the i — th rotor (Supplementary S4 and 
Movie MSI) and computed the spin-spin correlations. 

We find that there are three stages in the spin re- 
laxation process. In stage I, corresponding to the first 
~ 0.07 s, the rotors break their initial axial symmetry, 
Figure [2jt, and correlations decay rapidly with a char- 
acteristic Coulomb time scale t c ~ 0.02 s, Figure [2ji, 
which is the shortest time scale in the lattice relaxation, 

with t c ~ ^ a ^'° dominated by internal Coulomb inter- 
actions for the relaxation of a rotor interacting with two 
neighbors (Supplementary S4) in the absence of damping 
and external torques (inset of Figure HJi). Next, magnets 
of sub-lattices 1 and 2 (FigQJ,) organize in head to tail 
chains along the y direction, while those belonging to 
sub-lattice 3 still remain non-planar, Figure ^jp. In Stage 
II, once the sub-lattice 3 becomes planar, all the rods 
spin continuously leading to a plateau in the spin cor- 
relations (Fig[2}i); eventually the kinetic energy of the 
rotors has been dissipated sufficiently that the rotors os- 
cillate rather than spin. For our experimental parameters 
(FigQJi), the phase space trajectory changes from libra- 
tions to damped oscillations after 0.45 s (Supplementary 
Fig. S7); the rotors typically average about four full rota- 
tions before they switch to oscillations. Finally, in Stage 
III (Figl2j: and Fig|2ji) the rods oscillate without full ro- 
tations: when we fit the experimental dynamics at this 
state to a decaying exponential we find td ~ td, thus this 
stage is dominated by dissipative effects. 

To understand these different dynamical regimes, we 
performed molecular dynamics simulations of the massive 
undcrdamped rotors interacting through the full long- 
range internal Coulomb interactions between all the rods 
in the lattice using a Vcrlct algorithm (Supplementary 
S5, Fig.SlO and Movie SM2). In Figure EH, we see that 
the computed nearest neighbor spin correlations for the 
relaxation of the numerical lattice has the same three 
qualitative different regimes as in the experiments when 
the lattice relaxes from a polarized state to its spin ice 
manifold. Furthermore, the Coulomb and damping time 
scales for stage I and III as well as the plateau featuring 
stage II are in good agreement with experiments. The ob- 
served high frequency fluctuations in (S a Sp)(t) in both, 
experiments and simulations, are due to the Coulomb 



coupling between rods which rapidly reorient while they 
relax due to the fluctuations in the internal magnetic 
field. 

Having examined the dynamics of relaxation to the 
spin ice state, we now turn to the lattice response when 
a dipolc with charge \Q e \ at each pole and length L e , 
at a vertical distance h underneath the relaxed lattice 
is moved along one of the three sub-lattices at speed v 
(Supplementary Fig.S9). For an isolated rotor, the criti- 
cal torque that is required to destabilize the planar con- 
figuration is given by T c ~ 2aB c q, where B c is the ap- 
plied magnetic field; experiments on many rotors yielded 
an average B c ~ (2.4 ± 0.1) x 10~ 4 T. Equivalently, the 
threshold distance at which the external field will over- 
come both internal Coulomb interactions and static fric- 
tion is given by h* ~ y Q e Q a / 1 Tc- Dynamically, the in- 
ternal Coulomb interactions set a time scale for small 
out-of-plane oscillations of the rotors in the lattice, given 
by T p h ~ y A 2 I / [ioq 2 a ~ 0.01s for the experimental pa- 
rameters at hand. Thus, there are two dimensionlcss 
quantities that determine the response to the external 
perturbation: the ratio between phonemic and kinetic 
time scales vr p h/a and the ratio between internal and 
external magnetic forces, F mt / F ext = qh 2 /(Q e A 2 ). 

In Figure [3l we characterize the phase diagram of 
the dynamical response of the spin-ice lattice in terms 
of these dimensionlcss parameters. For h < h* , the 
lattice is disturbed only in a band of width D(h) ~ 
\J (h* 2 h) 2 / 3 — h 2 centered along the trajectory of the 
moving external dipole, based only on local interactions, 
static friction, and interactions with the external dipolc 
(Supplementary Fig.S9). For large v, T p h/rk > 1 so 
that the rotors have little time to respond and barely 
oscillate in an inertia-dominated regime. In the oppo- 
site limit, large amplitude oscillations and flips are ap- 
parent as there is enough time for the rotors to interact 
with the external dipole. Our results for these regimes 
show that the simulations (filled circles) and experiments 
(filled squares) agree. The solid line defines a threshold 
of the RMS fluctuations for the oscillations of all the rods 
^Li = 0-5 (Supplementary S4) separating the regimes. To 
understand this, we resort to a simple single rod approx- 
imation where the impulsive response of a rotor due to 
a long dipole located at a distance d(t) = y^h 2 + (vt) 2 , 
balances the change in its angular momentum yielding 
h ~ h*(Q e qa/I) 3 / 2 /v 3 , consistent with the observations 
when vTph/a 3> 1. Varying inertia from Iq to 4 Jo using 
our simulations we confirmed that as I grows, the bound- 
ary between interaction and inertial regime shift to the 
left; the inertial regime is reached for smaller values of 
vT ph /a and qh 2 /(Q e A 2 ). When h > h* , the Coulomb 
force due to the external field is not able to overcome the 
combined effects of static friction and internal Coulomb 
interactions, and the lattice falls into a friction dominated 
one in which oscillations are not apparent. 

Our spin ice phase emerges in a system of damped 
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macroscopic rotors, purely driven by interactions in a 
classical mechanical setting that differs from those found 
in its micro and mesoscopic relatives. Using a minimal 
model we can capture the dynamical evolution of the col- 
lection of rotors in the lattice observed in our experiments 
and reproduce the three main stages of lattice relaxation 
from a polarized state: explosive behavior lasting t c , dis- 
sipative librations and damped oscillations. The advan- 
tages of studying this macroscopic realization beyond the 
present work include the fact that (i) the interactions can 
be tuned through changes in the diameter of magnets 
or distance or orientation between them (Supplementary 
Fig. S4), (ii) incrtial and dissipative effects can be stud- 
ied by controlling the friction coefficient at the hinges 
as well as the mass of the rods, (iii) the effect of vacan- 
cies or random dilution can be examined by removing ro- 
tors from the lattice (iv) the lattice relaxation dynamics 
can be directly visualized at single particle level and (v) 
the system can be easily generalized to three dimensions 
(3-D) by stacking plates with hinged rotors along the 
z direction. Indeed a minimal 3-D realization is shown 
in Supplementary Fig. S12: a tetrahcdral configuration 
like the one found in the Pyrochlore lattice was created 
placing three acrylic plates one on the top of the other, 
the bottom and top plates contain three rotors defin- 
ing an equilateral triangle and the middle plate contains 
one rotor located equidistant from the others. The ease 
of fabrication, manipulation and measurement and the 
study of a variety of soft modes in artificial lattices in 
a system that is nearly five orders of magnitude larger 
and slower than its mesoscopic counterpart suggests that 
there is a new class of phenomena waiting to be explored 



in macroscopic frustrated systems. 
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FIG. 2: Lattice dynamics characterized by nearest neighbor spin correlations, {S a Sp). a, Stage I: once B z is turned 
off, the rotors originally pointing along z break their axial symetry. b, with the image showing the end of stage I and the 
onset of stage II when rods rotate respect to their center of mass yielding a plateau in (S a Sp). c, with a snapshot of the rods 
oscillating in stage III. d, In red experimental data obtained via image processing, in blue molecular dynamics simulation results 
from the numerical solution of equation (S4) where the full coulomb contributions from all neighbors is taken into account. At 
t=0 all S poles point along z. Stage I is dominated by Coulomb interactions between rods and characterized by the Coulomb 
time scale t c . In Stage II, all rotors spin until dissipation damps out the spin in favor of oscillations, leading to Stage III where 
they exhibit damped oscillations. After relaxation the rods lie in the x-y plane in a honeycomb spin ice magnetic configuration, 
with its characteristic nearest neighbor spin correlations {SaSp) = 1/3 (solid line). Inset: Experimentally measured value of 
(SaSp) during the initial explosive evolution (red) compared with cos(a) where a is the solution of Supplementary equation 
(S2) for one rotor interacting with two neighbors, in the absence of damping and external torques (green). 



6 




FIG. 3: Phase diagram of the lattice dynamical response to an external perturbation. The horizontal axis shows 
the dimensionless ratio of the kinetic and phononic time scales with v the speed of an external dipole, while the vertical axis 
shows the dimensionless ratio of the internal and the external magnetic forces due to an external dipole of strength Q located 
at distance h from the lattice (see text for details). Experimental and numerical data shown in squares and circles respectively, 
and colors define the nature of the lattice dynamical response to the external perturbation. We see that the dynamics may 
be broken up into a frictionally dominated, interaction dominated or inertially dominated regime as a function of the relative 
magnitude and rate of external forcing. 
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SI. Coulomb like interaction between two magnetic rods 



In order to compute the magnetic charge parameter q used in our simulations, we mea- 
sured the force-distance relation for several bar magnets using an Instron 5544 system to 
estimate the charge (pole strength) at the end of each bar in the dumbbell model. Two 
magnetic rods of equal dimensions were aligned with respect to each other at a distance z 
along the vertical axis and with opposite poles facing each other. Then, one of the rods was 
pulled apart and the interaction force measured, with a typical outcome shown in Figure 
1ST] where the experimental data for two rods in the configuration shown in the inset of 
Figure [HH is shown in red. The continuous black line is the fit obtained from a dumbbell 
model, where interactions between magnetic poles follow a Coulomb law, and the only free 
parameter is the magnetic charge at the tips of the rod: 

F ( Z ) = ^7 (~2 - ( z + 2a) 2 + ( 2 + 4a ) 2 ) ( S1 ) 

From the fit F(z) we found q = 2.03 ± 0.08 A m, in good agreement with previous 
estimates ]]. For a rod of radius r = d/2, q = M s iid 2 /4 [lj], using this result we computed 
the saturation mag netization as M s = (1.21 ± 0.03) x 10 6 A m" 1 which is in agreement 
with the available data of the magnetization of Neodymium rods, validating the dumbbell 
approximation when the rods are separated by a distance z > d, as previously pointed out 

mi 



S2. Moment of Inertia, Static friction and Damping coefficient. 

The moment of inertia of magnetic rods around their local rotation axis is I = ML 2 /12, 
where M is the rod mass (assuming uniform density) and L its length. The extra moment 
of inertia given by the plastic holder of mass m, is negligible as it is short and light m ~ 
0.05 x 10 -3 Kg. The mass of the magnetic rods used was measured individually for 50 
random rods which yielded: M = (0.278 ± 0.003) x 10 -3 Kg. From this data we computed 
Io = 8.41 x 10~ 9 Kg m 2 . This value is used as input in the numerical simulations. 

We quantified the static friction on each rotor, by placing a single rod at the center of a 
Helmholtz coil, as illustrated in Figure [S2] and measured the critical field at which the rod 
deviates from its initial position to find B c = (2.4 ±0.1) x 10~ 4 T. Using B c we obtained the 
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critical torque for which the rod starts to move and found that T c = 2aqB c ~ 0.93 x N 
m. 

To compute the damping coefficient 77, we isolated a single rotor and impulsively applied 
a torque to it, and then recorded its relaxation dynamics using a Phantom V7.3 high speed 
camera with frame rates between 1000 fps (frames per second) up to 4000 fps. Using 
standard imaging techniques we extracted the evolution of a(t), which corresponds to a 
damped dynamics in absence of external forcing. The damping is computed directly by 
fitting it to the solution a ~ exp(— £/td), (Fig. IS3l) and thus we estimated the damping 
time of a single rod to be r D = 0.83 ± 0.18 s. 



S3. Single rod Dynamics and Triad 

A magnetic bar of length L = 2a and diameter d <C 2a, behaves like two magnetic 
monopoles of equal strength but opposite sign located at each end of the bar {2), Oj. The 
magnetic moment of each magnet is then given by |m| = 2aq, where 2ar is the vector 
separating the two poles and m points from S to N, so that poles of two distinct rods 



interact through Coulomb's lawj^, (J. The dynamics of a single rotor can be understood by 
examining its response to a uniform magnetic field along the z direction, with an equation 
of motion given by: 

d 2 a da 

^ = Te " ^Tt (S2) 

where / is the rotational moment of inertia of the rod (Supplementary S2), 77 is the damping 
coefficient of the hinge, and T e = 2aqBo sin (a) is the torque due to the action of an external 
magnetic field on the localized charges at the end of the rod. This simplified system has two 
natural time-scales: an inertial time tb = 27r(//2a|g-Bo|) 1 ^ 2 ~ 0.003 / 1 1 1//2 an d a frictional 
time td = I /f] ~ Is. Since t# = td for a critical field Bq ~ 10~ 6 T, the dynamics of a single 
rod is underdamped in all our experiments. 

The minimal model where Coulomb interactions produces frustration in our system arises 
in the unit cell of three magnetic rods of length L at 120°, that is frustrated in the plane. 
When A/L < 0.2 the rods lie in the x — y plane in a spin ice configuration as shown in 
Figure [SU as A/L is increased, the out of plane magnetization of the triad M z increases 
(Fig. |S4]) . while when we reverse this operation, hysteresis is observed, a consequence of 
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static friction in the system (see Fig. ISllb for M z versus A/L in the lattice). We have used 
a single triad to examine the role of geometrical disorder by measuring M z when one of the 
hinges was rotated in the plane by 59 (Figure ED inset). We observe that the system remains 
in the x — y plane for small 59 ; however, when 59 ~ 35° M z starts to increase rapidly until 
59 ~ 50° when one magnet becomes perpendicular to the plane. 

The relaxation dynamics of the lattice made of unit cells is strongly dependent on the 
strength of the interactions and thus on the lattice parameters. For example, the damping 
time Tp iad of one of the rotors belonging to a triad grows with A (Fig. IS5l) . When the 
A/ 2a approaches the lattice parameters, the dynamics of one rod in the triad approaches 
that of a rod in the lattice. Indeed when A/2a < 0.2, we found r^ iad ~ 0.04 s, close 
to the relaxation time of a single rotor in the lattice, t re \ ~ 0.034 s, which was found by 
computing the experimental single-particle autocorrelation function C(t) = ((mj(t)mj(O)} — 
(rrii) 2 ) j \{m 2 ) — (rrii) 2 ) averaged over all rotors, and fitting its first 0.08 seconds of evolution 
to an exponential decay as shown in Figure IS6l 

S4. Experimental Lattice 

All 352 rods were made out of Neodymium (NdFeB plated with NiCuNi), with saturation 
magnetization M s = 1.2 x 10 6 A m _1 , length L = 2a — 1.9 x 10~ 2 m, diameter d = 1.5 x 10~ 3 
m and mass M = 0.28 x 10~ 3 Kg. They were hinged at the plane that is equidistant from 
their N and S magnetic poles, such that their axis of rotation crosses their center of mass. 
The hinges were made out of Acrylic-based photopolymer FullCure720 (Transparent) using 
a 3D printer Connex500 from ObJet Geometries. The hinges supporting the rods were 
introduced at the holes of an acrylic plate of size 0.61 x 0.61 m, defining the sites of a 
kagome lattice with lattice constant I = \^3(a + A) where A = 4.675 x 10~ 3 m. There was 
a small amount of geometrical disorder in the azimuthal orientation of the rotors, 9, due to 
lattice imperfections 59 max ~ 2°. We measured the azimuthal deviations of the rotors in the 
lattice when they were in the planar (x — y) configuration from lattice pictures (see below) 
to find that they follow a Gaussian distribution with mean 59 = 1.2° 

The number of rods, and their positions in the lattice were unchanged between exper- 
iments otherwise mentioned. To examine the configurations of the system, we took high- 
resolution pictures out of ten different experimental realizations, and used standard imaging 
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techniques to obtain 2D maps with the positions of each magnetic pole (South poles col- 
ored black). Each realization was accomplished after perturbing the lattice with an external 
dipole in order to produce a different initial and final magnetic state. In every occasion all 
vertices satisfied the ice rule. Away from the edges of the system, the spin ice rule was very 
robust, though small out of plane deviations larger than 5a, were sometimes observed at the 
edges of the lattice, a feature that becomes increasingly irrelevant for large systems. The 
averaged value of the total normalized magnetization of the sample, along x and y directions 
typically achieved (m x ) ~ 0.12 and {m y ) ~ —0.182 respectively. 

Correlations 

In order to obtain the evolution of the nearest neighbor magnetic correlations in time, 
rotors were polarized using a Solenoid of 850 Turns with length 0.1 m and radius r = 0.4 
m made out of copper wire of diameter d = 0.51 mm. The lattice was located inside the 
solenoid at its medial plane so that the radial component of the field was zero. The coil 
was connected in series with an ammeter and a regulated dc power supply serving as a 
current source. We used currents up to i = 3.5 A to create a magnetic field strong enough 
to polarize our sample. At i = 1.5 A, (B z ~ 3.2 x 10~ 3 T), all the N poles of the rods 
were pointing along the z direction. Next we turned the field off and waited until the rotors 
relaxed, simultaneously recording the evolution of the rods with a high speed Phantom V9.0 
camera. We then analyzed these data using standard imaging techniques and extracted the 
full time trajectory, ai(t), of each rod i from movies like the one shown in Supplementary 
SMI. 

The three dimensional evolution of the rods is characterized by the vector 

fhi (t) = (cos 6i sin (t) , sin 9i sin (t) , cos (t ) ) . 

thus the evolution of nearest neighbor rods correlations reads S a p(t)(rhi ■ rhi + i)(t) = 
(1/n) J^i^iit) ■ fni+i(t)/\mi(t) ■ m i+ i(t)\. A coarse grained charge at vertex k, at any time 
t, can be defined by Qkif) = J2i=i <li cos ot-iif)- First, second and third nearest neighbors 
magnetic correlation as well as nearest neighbor charge correlations after relaxation are 
summarized in Figure [S8l 
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Lattice dynamics 

The fastest time scale of lattice relaxation is dominated by Coulomb interactions at the 
very onset of motion where all rods are in a configuration with their S poles pointing along 
the z direction. During the first 0.07 seconds, two sublattices reach the plane before the 
remaining sublattice and organize in head to tail chains showing trend to ferromagnetic order. 
The experimental nearest-neighbor correlation decays exponentially, with a characteristic 
Coulomb time scale t c = 0.02 seconds. To understand this, we compute a(t) for a single rod 
using equation IS21 where torques are due to Coulomb interactions with neighbors oriented at 
120° relative to each other in the chains along the y direction. The numerical solution of this 
minimal model allows us to find that cos(a(i)) decays exponentially with characteristic time 
equal to t c . At leading order in a, it also gave us a good estimate for t c = ^nla/ /xo = 
0.02 s. 

Once all the rods reach the x — y plane, the next relaxation stage involves the rotation of 
the rods around their center of mass. The amount of time rotations last in the system can 
be found by numerically solving the nonlinear differential equation of a damped rotor which 
rotate due to the torque generated by Coulomb interactions with its four nearest neighbors 
located at 120° of each other in a mean field approximation: 

a" it) = Q 2 b sin(a(t)) - Q v a'{t) (S3) 

, with initial conditions a(0) = 0, a'(0) = V, where V is obtained from energy conservation 
at the onset of the rotor dynamics, fi& = \j2qa(B) rms , — r)/I ~ 1 seconds and (B) rms is 
the internal magnetic field due to the Coulomb interaction with its neighbors averaged over 
a cycle. For our experimental parameters, the phase space trajectory changes from open 
orbits into a dissipative attractor after 0.45 seconds: during this time the rotors average 
about 4 full rotations before they begin to oscillate as shown in Figure [S7l As expected, 
this simplified model does not reproduce all details of the collective dynamics. However, 
it captures well the typical time scales for full librations and oscillatory behavior expected 
before the collective ground state is reached. To resolve the effect of collective dynamics 
we solved the equations of damped motion for each rotor Coulomb interacting with all the 
other rotors in the system using molecular dynamics simulations as described in section S5. 
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Phase diagram 



To characterize the response of the system to external perturbations, we built a cart that 
can carry dipoles of length L e = 0.135 m, cross section radius r e = 0.00635, 0.00318, 0.00159 
m and charges Q e = 64g, 16q and 8q respectively. This cart moves underneath the lattice 
with speed v between 1 — 7 m/s, set by the initial impulse provided to the device (inset 
Fig. IS9a). The second experimental variable was the distance between the closest point 
of the dipole to the lattice, h, which varied between 0.42 to 0.05 m. In each realization 
we start with a different magnetic configuration where rotors lie at the x — y plane at not 
apparent magnetic order. The magnetic configuration is a result of agitating the relaxed 
lattice moving an external dipole randomly underneath the lattice. Then from the middle 
of the left boundary of the sample, a dipole with magnetic charge Q e and length L e and of 
strength B ext ~ Q e /h 2 , crossed the sample at a speed v from one lattice edge to the other. 
In order to distinguish between inertial and interaction regime, we quantified the response of 
the system computing the root mean square fluctuation (RMSF) in a, about the rod's x — y 
equilibrium position. This Lindemann like index 5n = ~ m ■ ™' ^ ~ > quantifies 
the angular fluctuations in the context of the parameters of the system, where ()t means 
the temporal average. The threshold value of the Lindermann parameters Sfy = 0.5. When 
8u < ^Lii the system is in an inertial regime. The boundary between interaction and inertial 
regimes can be obtained using a single particle approximation by equating the impulse a 
rotor feels when the external dipole is moving at speed v, at a distance d(t) = \Jh 2 + (ft) 2 , 
from it, to the change in the angular momentum of the rotor: 5L = f(F(t) x a)dt, where 
F = fioqQ e /{h 2 + {vt) 2 ). Integrating over time in the right side between —/3/v and /3/v, 
where (3 = J (h* 2 h) 2 / 3 — h 2 , yields the change in the angular momentum of the rod due 
to its interaction with the upper charge of the dipole: 5L ~ fi Q qQ e a/ (2whlv) (h* is the 
threshold height for which the Coulomb interaction between a rod and the external dipole 
overcomes the static friction, as explained in the paper). For an arbitrary size of oscillations, 
dfy, we can get an estimation of how h change with v by integrating over time at the left 
side of the previous equation to get : h ~ h* (fioQ e qa/ 1) 3 / 2 / v 3 , which gives good account 
of experimental results. The interplay between the ratio of internal to external magnetic 
forces and the speed of the dipole, determines in average the oscillations amplitude in a belt 
of the sample whose width D can be computed in the quasi-static approximation (Fig. IS9l) . 



7 



in terms of local interactions, static friction, and interactions with the top charge of the 
external dipole yielding D/h* ~ ^(h/h*) 2 / 3 - (h/h*) 2 , when h < h*. 



S5. Molecular Dynamics simulation 

Numerical results were obtained by direct numerical integration of the equations of motion 
for each rod interacting with all the others 351 rotors via Coulomb interactions in the 
dumbbell picture using a Verlet method with an integration time step At = 10~ 4 at zero 
temperature. In all simulations experimentally measured parameters for lattice constant, 
damping, inertia and charge of the rotors were used. Equilibration was accomplished after 
2 seconds of real time lattice evolution (2 x 10 4 time steps). To check the system size 
dependence on the calculated quantities, we compared our results with those for a system 
with n = 900 rods. As there was no significant difference in the results between the n = 900 
and n = 352 systems, we used the later number for all our simulations. 



T d 2 a® da® , 

~dtf~ = e ~ ' e: ^ ^ 

T (i) _ ( f cm rg 

e ~Utt )h 3 l^l 3 
fij is the vector joining charges qi and qj, and f^ m is the unit vector pointing from the center 
of mass (cm) of the rod where charge qi belongs. We find that the numerical lattice always 
fulfills the spin ice rules as shown in Figure IS101 Correlations after equilibration compared 
with the experimental values (Table of Fig. 158]) . show that nearest and second neighbor 
correlations compare well while third neighbors and charge correlations are smaller than 
in the experimental counterpart. We attribute this to the absence of static friction in the 
numerical model. To complement the experimental phase diagram we performed simulations 
using a point-like magnetic charge Q e = 64q, moving in a slightly larger parameter space that 
the one experimentally explored , E = (0.005,10) x (0.05,0.5). Our numerical constructed 
findings in the inertial and interaction regimes fall inside the experimental one, showing that 
the basics physics is captured well by the numerical model. 

Molecular dynamics simulations also permitted the examination of the out of plane lattice 
fluctuations versus variations in A/L. Figure ISllb shows the lattice RMS deviations out 
of the x — y plane averaged for all the n rods of the lattice M 2 rms = ^/Xa( cos «i 2 )/n 
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growing linearly for A < 0.3L. They show a behavior qualitatively similar to the total M z 
of three dumbbells interacting via Coulomb, in a star configuration, obtained from energy 
minimization, Figure ISllb for A < 0.4L. An extended numerical study will be published 
elsewhere. 



S6. Theory: Dipoles v/s Dumbbells 

Once all rotors were placed in the kagome lattice, we found that the honeycomb spin ice 
phase is very robust: any time a rotor located at the bulk was displaced with our finger tips, 
it came back to the x — y plane after a few oscillations. When we displaced it strongly so that 
a rotor flipped, locally new configurations in the spin ice manifold formed from sequential 
first-nearest-neighbor flips alone. Sequentially flipped second-nearest neighbor pairs were 
never observed, unless linked by a shared first nearest neighbor. Rotors being part of a 
closed magnetic loop were particularly stiff to external perturbations thus; single-particle 
spin-flipping suggests that flipping dynamics depends on local magnetic configurations. 

We examined two limit cases in order to understand these findings. In the first one, we 
let the magnets lie far apart (A/a > 0.4) so that a Hamiltonian with nearest neighbors 
dipolar interactions is a good approximation {2 — 4 1 . We found that the dipoles having a 
continuous U(l) symmetry, along their local axis prefer configurations in the x — y plane 
in the honeycomb spin ice manifold. The ground state energy is degenerate, E = —^-N t , 
where N t is the total number of triangles of the kagome lattice and D ~ 0(a/A) ~ 10 -5 J, 
sets the energy scale for interactions. The cost of raising a rod is SE = decreasing as 
A increases. This model also provides an easy estimation of the local energy configurations 
for a triad of rotors shown in Figure [Si] where the ones satisfying the spin ice rule are 
energetically favorable since head to tail magnetic configurations are energetically favorable 
(— 1.75D instead of 1.75D of tail to tail or head to head configurations) which explains also 
why magnets taking part of closed loops are particularly stable to perturbations. 

In the opposite limit, relevant for our experiments, rotors are brought close to each other 
such that A/a < 0.4, and thus a dumbbell model is a reasonable description The ground 
state of our system can be understood in terms of magnetic dumbbells with magnetic charges 
at the ends of each rod. For a single triad of dumbbell dipoles, each having a length 2a 
and a charge of magnitude ±g on each dipole head and tail respectively, the full Coulomb 
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interaction between magnetic charges contains 12 terms and an arbitrary constant that 
represents the self energy of this system. The minimal distance between the origin and 
the closest charge is A, therefore in the limit A — > the energy will be given only by the 
three divergent Coulomb terms due to the interactions between the three charges 1,2,3, 
that are closest to the origin. Thus, at leading order, the interaction energy is given by 
Ua ~ U12 + U23 + U31, where is clearly divergent. In the same spirit as in QED and QCD 
5|, we describe the basic physics by paying attention only to the most divergent quantities 
and can understand the optimal charge configuration using the fact that ~ fl'/A. Thus 
minimizing the Coulomb energy is equivalent to minimizing g. In the limit when A — y 0, a 
is the only length scale. Then the leading order term in g will only depend on the geometry, 
and any needed change in the sign of the charge product will appear as a it change in 
the polar rotation angle. Therefore in considering the lowest energy configuration of three 
interacting dumbbells that are at a distance A from each other, we find that: 

g(a 1 ,a 2 ,a 3 ) = 
(5 + cos «i cos «2 — sin ot\ (—3 + sin a 2 ) — 3 sin c^)^ 1 / 2 + 
(5 + cos ot\ cos 03 + sin a\ (3 + sin a 3 ) + 3 sin a 3 ) _1 / 2 + 
(5 — cos 0L2 cos «3 — sin «2 (3 + sin a 3 ) + 3 sin a^)" 1 / 2 

where a; are the rotation angles measured relative to the vertical axis (polar angle). The 
global minimum for g(a\, «2, 03) occurs for planar configurations corresponding to the spin 
ice case, i.e. (a±, 02,03) = ±(7r/2, it/2, — it/2) or permutations. Thus, the basic building 
blocks of our experiment follow spin ice rules as the shortest distance between magnetic 
poles goes to zero asymptotically: the initial U(l) symmetry for each rotor is reduced into 
a Z2 Ising like symmetry. 



Movies 



Caption Movie SMI 

Experimental lattice dynamics. In the initial state all the S poles of the lattice n 
rotors are pointing along z. At t = 0, the magnetic field is switched off and the lattice 
relaxes over a period of 2 seconds. This dynamical process is shown here along with the 
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evolution of the nearest neighbor magnetic correlations. In the final state all the magnets 
lie in the x — y plane in a honeycomb spin ice magnetic order. 

Caption Movie SM2 

Computational lattice dynamics. In the initial state all the S poles of the lattice 
n rotors are pointing along z. At t = 0, the magnetic field is switched off and the lattice 
relaxes over a period of 2 seconds. This dynamical process is shown here along with the 
evolution of the nearest neighbor magnetic correlations. In the final state all the magnets 
lie in the x — y plane in a honeycomb spin ice magnetic order. 
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FIG. SI: Force vs. distance between two Neodymium rods. In red are shown the 
experimental data with error bars of the attractive force between two magnetic rods, as a 
function of the distance between the two closest faces of the rods. Black line represents the 
best fit obtained by using the Coulomb law (Eq. flSlj) ). from where we obtained the 
magnitude of the magnetic charges q = 2.03 ± 0.08 A m. The inset illustrates a typical 

tensile experiment. 
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FIG. S2: Schematic of the setup used to measure static friction. A uniform 
magnetic field is applied to a single rod producing a torque in each of its poles. Once the 
rod depart from its equilibrium position we recorded the value of the magnetic field and 

computed T c . 
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FIG. S3: Single rod dynamics. In blue are the experimental data taken from standard 
imaging techniques and the dashed black is the best fit obtained using an exponential 
damping model, offset by 0.1 for clarity purposes. 
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FIG. S4: A and 6 dependence of the normalized magnetization for three rods 
in a 120° star configuration. M z grows with A /2a and when the operation is reversed, 
the system shows hysteresis. The inset shows M z growing when one of the rotors change 
its azhimutal orientation respect to the others by 89 
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FIG. S5: r^ iad versus 2a/ A. The relaxation time of a single rotor in a triad 
configuration with two nearest neighbors. T^ iad ~ t re [, when 2a/ A < 0.2. The data 
2a/ A = is the damping time of a free rod, td. 
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FIG. S6: Single particle autocorrelation, C(t). The (blue) points were obtained by 

computing the experimental single-particle autocorrelation function 
C(t) = ((mj(t)mj(O)) — (mi) 2 ) / \(m 2 ) — (mi) 2 ) averaged over all rotors. Inset: C(t) for the 
first 0.08 seconds of lattice relaxation (blue) compared with the function cos(exp(— t/t re i)) 
which is the the best fit for C(t) (dashed black curve). From the fit we obtained the 
relaxation time of a single rotor in the lattice t re \ = 0.034 sec. 
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FIG. S7: Lattice dynamic at stage II. The phase space trajectory corresponding to 
the solution of equation [S3] for a single rotor interacting with its four neighbors, after 0.45 

seconds have elapsed. 
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FIG. S8: Correlations after relaxation. Top: The honeycomb structure formed by 
connecting the spins of the kagome lattice. Each bar element represents a rod oriented 
along the axis. The Greek symbols label spins for correlation calculations. Bottom: First, 
second and third nearest neighbors magnetic correlation and nearest neighbors charge 
correlations for the experimental and numerical lattices. 
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FIG. S9: D/h* v/s h/h*. a, Experimental data showing the width D of the stripe of 
the sample which is excited due to two dipoles of charges Q e = 16q (blue) and 64g (red), in 

good agreement with the theoretical scaling (black dots) and the results from simulations 
(green) for Q e = 64g. The inset illustrates the experiment, where a dipole of length L e and 

charge Q e located at distance h from the lattice is exciting a stripe of magnetic rods, b, 
Superimposed images obtained by computirj.g the intensity difference between frames when 

the lattice is excited by an external dipole and a typical image from where we obtained 
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FIG. S10: Numerical lattice: Rods in a typical configuration fulfilling spin-ice rules 
after the numerical lattice has reached relaxation. 
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FIG. Sll: Numerical M z v/s A. a, RMS deviations out of the x — y plane averaged for 
all the n rods of the lattice M z rms v/s A /2a after 4 seconds of simulation have elapsed. 
As the distance between nearest poles increases, rotors are more prone to choose positions 
which depart from the x — y plane at equilibrium, b, Average magnetization along the z 

direction of a triad of rotors v/s A /2a. M z was obtained computing the Coulomb 
interaction between the six poles and minimizing numerically the energy en the domain 

(0,2tt) 3 . 
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FIG. S12: A minimal 3-D generalization of our system. A tetrahedral 
configuration like the one found in the Pyrochlore lattice was created placing three acrylic 
plates one on the top of the other, the bottom and top plates contain three rotors defining 
an equilateral triangle and the middle plate contains one rotor located equidistant from the 

others. 
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